/*****************************************************
*         Università degli Studi di Cagliari         *
*                                                    *
*  Algoritmi e Strutture Dati 2 - Progetto mid-term  *
*                                                    *
*  Riccardo Scateni    (riccardo@unica.it)           *
*  Stefano Marras      (stefano.marras@unica.it)     *
*  Marco Livesu        (marco.livesu@unica.it)       *
******************************************************/

#include "core.h"

void Core::open_file()
{

        QStringList files = QFileDialog::getOpenFileNames( NULL,
                                                       "Open mesh",
                                                       ".",
                                                       "3D Meshes (*.off *.obj *.ply);;OFF (*.off);;OBJ (*.obj);;PLY (*.ply)");


        QStringList::Iterator it = files.begin();


        while( it != files.end() )
        {

                if( !(*it).isNull() ) loadMesh(*it);

                ++it;
        }
}



void Core::loadMesh(QString filename)
{

        MyMesh *m = new MyMesh();

        m->request_face_normals();
        m->request_vertex_normals();
        m->request_vertex_colors();


        OpenMesh::IO::Options options;

        if ( OpenMesh::IO::read_mesh(*m, filename.toStdString(), options) )
        {
                QFileInfo info(filename);

                OpenMesh::MPropHandleT<QString> filename;
                m->add_property(filename, "m:filename");
                m->property(filename) = info.fileName();


                if ( !options.check( OpenMesh::IO::Options::FaceNormal ) )   m->update_face_normals();
                if ( !options.check( OpenMesh::IO::Options::VertexNormal ) ) m->update_vertex_normals();


                // bounding box
                MyMesh::ConstVertexIter vIt(m->vertices_begin());
                MyMesh::ConstVertexIter vEnd(m->vertices_end());

                OpenMesh::Vec3f bbMin, bbMax;

                bbMin = bbMax = OpenMesh::vector_cast<OpenMesh::Vec3f>(m->point(vIt));

                ++vIt;

                for ( ; vIt!=vEnd; ++vIt)
                {
                        bbMin.minimize(OpenMesh::vector_cast<OpenMesh::Vec3f>(m->point(vIt)));
                        bbMax.maximize(OpenMesh::vector_cast<OpenMesh::Vec3f>(m->point(vIt)));
                }

                OpenMesh::MPropHandleT<BBox> bbox;
                m->add_property(bbox, "m:bbox");
                m->property(bbox).min = bbMin;
                m->property(bbox).max = bbMax;

                generateSpherical2DTextureMapping(m);
                setRandomFaceColor(m);
                setRandomVertexColor(m);

                mesh.append(m);

                emit updateMeshList(&mesh);

                return;
        }
}



void Core::meanCurvature()
{

        for (int i=0; i < mesh.size(); ++i)
        {
                OpenMesh::VPropHandleT<MyMesh::Scalar> curv = request_meanCurvature(mesh[i]);
                mesh[i]->remove_property(curv);
        }
}



void Core::GaussianCurvature()
{

        for (int i=0; i < mesh.size(); ++i)
        {
                OpenMesh::VPropHandleT<MyMesh::Scalar> curv = request_GaussianCurvature(mesh[i]);
                mesh[i]->remove_property(curv);
        }
}




OpenMesh::EPropHandleT<MyMesh::Scalar> Core::request_CotagentWeights(MyMesh *m)
{

        OpenMesh::EPropHandleT<MyMesh::Scalar> cotan;
        m->add_property(cotan, "e:cotan");

        MyMesh::EdgeIter eIt(m->edges_begin()),
                         eEnd(m->edges_end());

        MyMesh::HalfedgeHandle he0, he1, he2;

        MyMesh::Point    p0, p1, p2, d0, d1;

        MyMesh::Scalar   c, weight, lb(-0.99999), ub(0.99999);

        for ( ; eIt != eEnd; ++eIt)
        {

                weight = 0.0;

                he0 = m->halfedge_handle(eIt, 0);
                p0  = m->point(m->to_vertex_handle(he0));

                he1 = m->halfedge_handle(eIt, 1);
                p1  = m->point(m->to_vertex_handle(he1));

                if ( !m->is_boundary(he0) )
                {

                        he2 = m->next_halfedge_handle(he0);
                        p2  = m->point(m->to_vertex_handle(he2));

                        d0 = (p0 - p2).normalize();
                        d1 = (p0 - p1).normalize();

                        c = OpenMesh::dot(d0, d1);

                        if (c < lb) c = lb;
                        else
                        if (c > ub) c = ub;

                        weight += 1.0 / tan(acos(c));
                }

                if ( !m->is_boundary(he1) )
                {

                        he2 = m->next_halfedge_handle(he1);
                        p2  = m->point(m->to_vertex_handle(he2));

                        d0 = (p0 - p2).normalize();
                        d1 = (p0 - p1).normalize();

                        c = OpenMesh::dot(d0, d1);

                        if (c < lb) c = lb;
                        else
                        if (c > ub) c = ub;

                        weight += 1.f / tan(acos(c));
                }

                m->property(cotan, eIt) = weight;
        }

        return cotan;
}


//  Calculus of the per vertex Voronoi area, as defined in
//  "Discrete differential-geometry operators for triangulated 2-manifolds"
//  Meyer, Desbrun, Schroder, Barr - Visualization and Mathematics III
//  Springer-Verlag, Heidelberg (2003)
//
//  A_Voronoi(v_i) = 1/8 * \sum_{j \in N(v_i)} { (cot \alpha_{ij} + cot \betha_{ij}) * \Vert v_i - v_j \Vert }
//
OpenMesh::VPropHandleT<MyMesh::Scalar> Core::request_VoronoiArea(MyMesh *m)
{

        OpenMesh::VPropHandleT<MyMesh::Scalar> area;
        m->add_property(area, "v:area");

        MyMesh::VertexIter vIt(m->vertices_begin()),
                           vEnd(m->vertices_end());

        MyMesh::FaceIter   fIt(m->faces_begin()),
                           fEnd(m->faces_end());

        MyMesh::FaceVertexIter fvIt;

        MyMesh::VertexOHalfedgeIter voheIt;

        MyMesh::VertexHandle v0, v1, v2;

        MyMesh::Scalar a;

        OpenMesh::EPropHandleT<MyMesh::Scalar> cotan = request_CotagentWeights(m);


//        for ( ; vIt != vEnd; ++vIt)
//        {
//                MyMesh::Scalar VoronoiArea = 0.0;

//                if ( !m->is_boundary(vIt) )
//                {
//                        MyMesh::Point pi = m->point(vIt);

//                        voheIt = m->voh_iter(vIt);

//                        do
//                        {
//                                MyMesh::Point pj = m->point(m->to_vertex_handle(voheIt.handle()));

//                                MyMesh::EdgeHandle e = m->edge_handle(voheIt.handle());

//                                VoronoiArea += (pi - pj).norm() * m->property(cotan, e);
//                        }
//                        while ( ++voheIt );
//                }

//                VoronoiArea /= 8.0;

//                m->property(area, vIt) = VoronoiArea;
//        }

        for ( ; vIt != vEnd; ++vIt)
        {

                m->property(area, vIt) = 0.0;
        }


        for ( ; fIt != fEnd; ++fIt)
        {

                fvIt = m->fv_iter(fIt);

                v0 = fvIt.handle();  ++fvIt;
                v1 = fvIt.handle();  ++fvIt;
                v2 = fvIt.handle();

                MyMesh::Point p0 = m->point(v0);
                MyMesh::Point p1 = m->point(v1);
                MyMesh::Point p2 = m->point(v2);

                a = 0.5 * OpenMesh::cross(p1 - p0, p2 - p0).norm();

                m->property(area, v0) = a / 3.0;
                m->property(area, v1) = a / 3.0;
                m->property(area, v2) = a / 3.0;
        }

        return area;
}



OpenMesh::VPropHandleT<MyMesh::Scalar> Core::request_GaussianCurvature(MyMesh *m)
{
        MyMesh::VertexIter vIt(m->vertices_begin()),
                           vEnd(m->vertices_end());

        OpenMesh::VPropHandleT<MyMesh::Scalar> GaussianCurv;
        m->add_property(GaussianCurv, "v:GaussianCurv");

        MyMesh::VertexVertexIter vvIt;

        OpenMesh::VPropHandleT<MyMesh::Scalar> area  = request_VoronoiArea(m);

        for ( ; vIt != vEnd; ++vIt)
        {

                MyMesh::Scalar sum = 0.0;

                if ( !m->is_boundary(vIt) )
                {

                        MyMesh::Point p0 = m->point(vIt);

                        vvIt = m->vv_iter(vIt);
                        MyMesh::Point p1 = m->point(vvIt.handle());

                        ++vvIt;

                        do
                        {

                            MyMesh::Point p2 = m->point(vvIt.handle());

                            MyMesh::Scalar angle = OpenMesh::dot(p1 - p0, p2 - p0);

                            angle /= (p1 - p0).norm() * (p2 - p0).norm();

                            p1 = p2;

                            sum += angle;
                        }
                        while ( ++vvIt );

                        sum = (2 * M_PI - sum) * (1.0 / m->property(area, vIt));
                }

                m->property(GaussianCurv, vIt) = sum;
        }



        // Determine min and max curvature and remove outliers
        //
        m->request_vertex_texcoords1D();
        {
                MyMesh::VertexIter vIt(m->vertices_begin()),
                                   vEnd(m->vertices_end());


                QVector<MyMesh::Scalar> values;

                for ( ; vIt != vEnd; ++vIt)
                {

                        values.append( m->property(GaussianCurv, vIt));
                }

                qSort(values.begin(), values.end());

                unsigned int n = values.size();
                unsigned int i = n / 20;
                MyMesh::Scalar min = values[i];
                MyMesh::Scalar max = values[n-1-i];

                for ( vIt = m->vertices_begin(); vIt != vEnd; ++vIt)
                {

                        m->set_texcoord1D(vIt, MyMesh::TexCoord1D((m->property(GaussianCurv, vIt) - min) / (max - min)) );
                }
        }


        // clean-up properties
        m->remove_property(area);

        return GaussianCurv;
}



OpenMesh::VPropHandleT<MyMesh::Scalar> Core::request_meanCurvature(MyMesh *m)
{
        MyMesh::VertexIter vIt(m->vertices_begin()),
                           vEnd(m->vertices_end());

        MyMesh::VertexOHalfedgeIter vheIt;

        MyMesh::VertexHandle v;

        MyMesh::EdgeHandle e;

        MyMesh::Point laplace;


        OpenMesh::VPropHandleT<MyMesh::Scalar> meanCurv;
        m->add_property(meanCurv, "v:meanCurv");


        OpenMesh::EPropHandleT<MyMesh::Scalar> cotan = request_CotagentWeights(m);
        OpenMesh::VPropHandleT<MyMesh::Scalar> area  = request_VoronoiArea(m);


        for ( ; vIt != vEnd; ++vIt)
        {

                laplace = MyMesh::Point(0.0, 0.0, 0.0);

                if ( !m->is_boundary(vIt) )
                {

                        vheIt = m->voh_iter(vIt);

                        do
                        {
                                v = m->to_vertex_handle(vheIt.handle());
                                e = m->edge_handle(vheIt.handle());

                                laplace += m->property(cotan, e) * (m->point(v) - m->point(vIt));
                        }
                        while ( ++vheIt );

                        laplace /= 2.0 * m->property(area, vIt);
                }

                m->property(meanCurv, vIt) = 0.5 * laplace.norm();
        }

        // Determine min and max curvature and remove outliers
        //
        m->request_vertex_texcoords1D();
        {

                MyMesh::VertexIter vIt(m->vertices_begin()),
                                   vEnd(m->vertices_end());


                QVector<MyMesh::Scalar> values;

                for ( ; vIt != vEnd; ++vIt)
                {

                        values.append( m->property(meanCurv, vIt));
                }

                qSort(values.begin(), values.end());

                unsigned int n = values.size();
                unsigned int i = n / 20;
                MyMesh::Scalar min = values[i];
                MyMesh::Scalar max = values[n-1-i];

                for ( vIt = m->vertices_begin(); vIt != vEnd; ++vIt)
                {

                        m->set_texcoord1D(vIt, MyMesh::TexCoord1D((m->property(meanCurv, vIt) - min) / (max - min)) );
                }
        }


        // clean-up properties
        m->remove_property(area);
        m->remove_property(cotan);

        return meanCurv;
}


void Core::setRandomVertexColor(MyMesh *m)
{
        MyMesh::VertexIter vIt(m->vertices_begin()),
                           vEnd(m->vertices_end());


        for (; vIt != vEnd; ++vIt)
        {

                uchar r = rand();
                uchar g = rand();
                uchar b = rand();

                m->set_color(vIt, MyMesh::Color(r,g,b));
        }
}


void Core::setRandomFaceColor(MyMesh *m)
{
        MyMesh::FaceIter fIt(m->faces_begin()),
                         fEnd(m->faces_end());


        for (; fIt != fEnd; ++fIt)
        {

                uchar r = rand();
                uchar g = rand();
                uchar b = rand();

                m->set_color(fIt, MyMesh::Color(r,g,b));
        }
}


void Core::setRandom1DTextCoord(MyMesh *m)
{
        m->request_vertex_texcoords1D();

        MyMesh::VertexIter vIt(m->vertices_begin()),
                           vEnd(m->vertices_end());


        srand(time(NULL));

        for (; vIt != vEnd; ++vIt)
        {

                float coord = rand()%2560 / 2560.0;

                m->set_texcoord1D(vIt, coord);
        }
}


void Core::generateSpherical2DTextureMapping(MyMesh *m)
{
        m->request_vertex_texcoords2D();

        MyMesh::VertexIter vIt(m->vertices_begin()),
                           vEnd(m->vertices_end());

        for ( ; vIt != vEnd; ++vIt)
        {

                MyMesh::Normal normal = m->normal(vIt);

                MyMesh::TexCoord2D texcoord2D(asin(normal[0]) / M_PI + 0.5, asin(normal[1]) / M_PI + 0.5);

                m->set_texcoord2D(vIt, texcoord2D);
        }
}
